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Abstract 

The Monte Carlo with Absorbing Markov Chains (MCAMC) method is introduced. This method is a generalization 
of the rejection-free method known as the n-fold way. The MCAMC algorithm is applied to the study of the 
very low-temperature properties of the lifetime of the metastable state of Ising ferromagnets. This is done both 
for square-lattice and cubic-lattice nearest-neighbor models. Comparison is made with exact low-temperature 
predictions, in particular the low-temperature predictions that the metastable lifetime is discontinuous at particular 
values of the field. This discontinuity for the square lattice is not seen in finite-temperatures studies. For the cubic 
lattice, it is shown that these 'exact predictions' are incorrect near the fields where there are discontinuities. The 
low-temperature formula must be modified and the corrected low-temperature predictions are not discontinuous 
in the energy of the nucleating droplet. 

Key words: Advanced algorithms, Dynamic Monte Carlo, lifetimes, metastable, prefactor 
PACS: 05.10.-a, 05.10.Ln 



1. Introduction 

One of the most difficult problems facing sim- 
ulations in science and mathematics is to be able 
to simulate time and length scales comparable to 
those found in nature, those required for engineer- 
ing applications, and those accessible to experi- 
ments. In this paper a method to obtain extremely 
long simulations for discrete systems is reviewed. 
This method, which merges the ideas of absorb- 
ing Markov Chains with those of kinetic Monte 
Carlo simulations, is called the MCAMC algorithm 
(Monte Carlo with Absorbing Markov Chains) [1- 
3]. It expands on the ideas put forward by Bortz, 
Kalos, and Lebowitz [4] in a method called the n- 



fold way. The MCAMC algorithm with a single 
transient state (which corresponds to the current 
state of the configuration) is the discrete-time ver- 
sion of the n-fold way algorithm [5] . 

As an application requiring long-time simula- 
tions, consider the problem of thermal reversal 
of nanoscale ferromagnetic grains. For highly 
anisotropic materials each region of the ferromag- 
net may be considered to have two states corre- 
sponding to two discrete directions. This gives an 
Ising model as an approximate Hamiltonian for 
ferromagnetic nanoparticles [6]. Starting from a 
quantum Hamiltonian of spin ~ particles inter- 
acting with a heat bath, it is possible to derive in 
certain limits a physical dynamic for the model 
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[7,8]. This dynamic has two parts: first a spin is 
chosen uniformly from all the spins, then the spin 
is flipped with some probability pm P that depends 
on the temperature T and the applied magnetic 
field H. One algorithmic step is called a Monte 
Carlo step (mcs) and corresponds to one spin-flip 
attempt. For a system with iV Ising spins, the sim- 
ulation time is often given in Monte Carlo steps 
per spin (MCSS), which corresponds to N mcs. 

With coupling to a fermionic heat bath [7], the 
probability pfn p is the Glauber transition proba- 
bility [9]. In this case p flip = e -P E ™-~ /[ e -P E ™ + 
e -/3£ old j where (3 = T^ 1 (with Bolztmann's con- 
stant set to one) , -E n ew is the energy of the config- 
uration with the Ising spin flipped and E \d is the 
energy of the current spin configuration of the sys- 
tem. A different probability (which also satisfies 
detailed balance) for coupling to a <i-dimcnsional 
bath of phonons has been recently derived [8] . 

The reason long-time simulations are required is 
now readily apparent. The underlying algorithmic 
step corresponds to an inverse phonon frequency, 
about 10~ 13 seconds. The time scale for simula- 
tion must be on the order of years to decades for 
simulations applicable to magnetic recording - 
basically simulations at least as long as data in- 
tegrity for a written bit of information. For paleo- 
magnetism simulations, time periods of millions of 
years must be achievable. The typical clock speed 
of a computer is 1CP 9 seconds. Consequently, 
faster-than-real-time simulations are required 
for realistic feasible simulations. Note that these 
faster-than-real-time algorithms cannot change 
the dynamic which has been derived starting from 
the quantum Hamiltonian. Hence advanced algo- 
rithms that are used in static simulations, such as 
reweighting techniques, the Swendsen-Wang algo- 
rithm, multicanonical methods, or simulated tem- 
pering [3] cannot be used — they would change 
the underlying dynamic. Rather, the physical dy- 
namic described above must be implemented on 
the computer in a more intelligent fashion than a 
brute-force method. 

This paper describes only the MCAMC al- 
gorithm and results obtained from it for low- 
temperature Ising simulations. Two other recent 
advances using rejection-free methods should be 
mentioned [3]. One is to use massively parallel 



computers in such simulations [3,10-12]. Another 
is to generalize the rejection-free algorithms to 
continuous systems, such as Heisenberg spin sys- 
tems [13]. 

Sec. 2 describes the MCAMC algorithm for 
the simple cubic (sc) nearest-neighbor (nn) Ising 
model. The generalization to other discrete sys- 
tems is straightforward. In Sec. 3 low-temperature 
predictions for the Ising model are reviewed. Sec- 
tions 4 and 5 present results for low-temperature 
metastable lifetimes in d=2 and d=3. Sec. 5 con- 
tains conclusions. 



2. The MCAMC Algorithm 

We simulate the nn Ising model with Hamilto- 
nian 7i = —Jj2(i j) a i (J j — H J2 i <Ji with J>0 the 
ferromagnetic coupling constant, the first sum run- 
ning over all nn pairs, and the second sum running 
over all N spins. We consider only the case of pe- 
riodic boundary conditions. The simulation starts 
with all spins <7j = 1 (all spins up), and a field H < 
(external field down) is applied. The magnetiza- 
tion is M = J2i °i ■ We measure the time r to go 
from the starting state to a state with M = 0. This 
simulation is repeated many times (for this paper 
1000 times) for the same T and H values to obtain 
the average lifetime (t) of the metastable state. 

For the sc lattice for every configuration each 
spin can be classified to be in one of 14 possible 
classes. (For the square-lattice Ising model there 
are 10 classes.) The spin is either up or down, and 
we label the first 7 classes with spin up and classes 8 
through 14 with spin down. The other factor deter- 
mining the class of a spin is the number of nn spins 
that are up, which can be any integer from zero to 
6. Let class 1 have 6 nn spins up, class 2 5 nn spins 
up, etc. Class 8 has 6 nn spins up, • • •, class 14 has 
nn spins up. Let rn be the number of spin in class 
i. Then N = J2iti n i- Let Pi be the probability of 
flipping a spin in class i, given that the spin was 
chosen during the first part of the dynamic algo- 
rithm. Then the probability of flipping any of the 
spins in class i in one mcs is niPi/N since rn/N is 
the probability of choosing a spin from class i. 

In order to exit from the current spin configura- 
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tion, one of the spins in one of the 14 classes must 
be flipped. Define Qi — jj S}=i n jPj f° r eacn of 
the 14 spin classes, and Qo = 0. Let [-J denote the 
integer part of a number. The discrete n-fold way 
algorithm has three steps. First, the time m to exit 
from the current configuration (in unit of mcs) is 
determined by 



m 



ln(r) 



ln(l - Q u ) 



+ 1 



(1) 



with f a random number uniformly distributed be- 
tween zero and one. Then using a different uni- 
formly distributed random number f the spin class 
k is chosen such that it satisfies Qk~i < rQu < 
Qk- Finally, using a third uniformly distributed 
random number, one of the spins in class k 
is chosen, and this spin is flipped. Of course, the 
flipped spin changes its class, as do its 6 nn spins. 
This means that the numbers rij change after one 
algorithmic step. 

The first two steps of the discrete n-fold way 
algorithm correspond to having an absorbing 
Markov chain with 14 states in the recurrent ma- 
trix with elements given by the probability to exit 
from the current spin configuration by choosing 
and flipping a spin in class i, namely riiPi/N. The 
transient subspace of the absorbing Markov ma- 
trix is a scalar (s = 1) given by the matrix T s — i = 
N(l - Qu). The time in the s = 1 MCAMC simu- 
lation corresponds to the time required to leave the 
current spin configuration. At low temperatures 
this time may be extremely long (if Q14 is very 
small), and the algorithmic speed-up obtained in 
the simulation can be many orders of magnitude. 

The simulation starts with all spins up. In one al- 
gorithmic step the s = 1 MCAMC described above 
exits from this state with a probability p\ (since all 
spins are up, n\ — N for this state). For low tem- 
peratures and H < 2J the most probable thing 
that happens in the next algorithmic step is that 
the flipped spin is again chosen and the next config- 
uration is again all spins up. This occurs with prob- 
ability ratio ( jv ~ 7 )P 1+6 P 2 _ By increasing the size 
of the transient surjspace of the absorbing Markov 
chain, it is possible to exit the absorbing Markov 
chain in a configuration with more than one spin 
flipped. For example, the absorbing Markov matrix 
to exit from the state with all spins up into a state 



with two spins flipped has two states in the tran- 
sient subspace (s = 2) with the transient matrix 



N \N Pl N(l- Pl ) 
and the recurrent matrix 

_ 1 [QP2 (N-7) Pl 



N V 







(2) 



(3) 



withx2 = N—pg—6p2 — (N—7)pi. The ratio to exit 
to a state with two nn spins overturned to that to 
exit to a state with two non-nn spins overturned is 
6p 2 /{N — 7)pi. The vector representing the initial 
state of all spins up is if J = ( 1 ) , and the 
probability to be in each of the two transient states 
after one time step is if JT S=2 ■ The first transient 
state (lower right-hand corner of the matrix) is the 
state with all spins up, and the other state (upper 
left-hand corner of the matrix) corresponds to the 
N states that have one overturned spin. By adding 
more states to the transient subspace longer times 
to exit the transient subspace are possible. 

The time m (in mcs) to exit from the transient 
subspace is found from the solution of the equation 



if < r < 



lfjT n 



(4) 



where if is the vector of length s with all elements 
equal to unity. The probability that the system 
ends up in a particular absorbing state is given by 
the elements of the vector T^NR with the funda- 
mental matrix N = (I — T) _1 . 

To obtain very long lifetimes, a multiple preci- 
sion package was used [14]. To keep bookkeeping 
overhead to a minimum, for the Ising simulations 
below, only s = 1 MCAMC was used except when 
the spin configuration had all spins up, then higher 
s MCAMC (up to s = 5) was used. 



3. Low-temperature Metastable Ising 
Predictions 

At very low temperatures the kinetic Ising sim- 
ulations are influenced by the discreteness of the 
lattice. This allows for the exact calculation of the 
saddle point as well as the most probable route to 
the saddle point [15]. For the square lattice [15] 
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with £ 2 = [2J/\H\\ + 1 the average lifetime (in 
units of mcs) at low temperature is given by (r) = 
A 2 exp(T 2 /T) with r 2 = 8J£ 2 - 2\H\(t% - £ 2 + 1). 
The prefactor, A 2 , was determined to be 5/4 for 
l 2 = 1 and 3/8 for i 2 = 2 [2]. Recently, the pref- 
actor has been found to be A 2 = 3/ [8 (£ 2 — 1)] for 
all \H\ < 2 J [16]. Hence the nucleating droplet for 
the square lattice is a rectangle of overturned spins 
of size l 2 x (£ 2 — 1) with one overturned spin on 
a long interface. At low enough temperature, this 
result for (t) should hold everywhere, except when 
2 J/\H\ is an integer. Note that the prefactor is dis- 
continuous for the special values where 2J/\H\ is 
an integer. 

For the sc lattice, the average lifetime is given by 
<r> = A 3 exp(r 3 /T) with £ 3 = [4J/\H\\ + 1 and' 

r 3 = Y1J£\ - 8J£ a - 2\H\£ 2 3 (£ 3 -1)+T 2 (5) 

[16-18]. This result should hold when 4J/\H\ 
is not an integer and when £3 > 2. The pref- 
actor has recently been found to be A3 = 
[16(4 - h + 1)(£ 2 - 1)] _1 [16] for £ 3 > 2. Hence 
the nucleating droplet is a cube with length £3 
with one layer removed and the square-lattice nu- 
cleating droplet placed on this layer. Note that at 
the special values where AJ/\H\ is an integer T 3 is 
discontinuous, and in fact decreases by 4J. 



4. Square-lattice Results 



The MCAMC values obtained for the square- 
lattice lifetime at \H\ = 3/4 (so £ 2 = 3) are shown 
in Fig. 1. Note that the average lifetime is very 
large. For L = 20 if one Monte Carlo step per 
spin (MCSS) corresponded to about 10~ 13 sec then 
the largest time corresponds to about one year. 
The low temperature prediction with prefactor one 
(dashed line) does not describe the simulation data 
very well. However, including the low temperature 
prefactor [16] of 3/16 for this field value (heavy 
solid line) fits the data very nicely. 

Figure 2 presents results for Tln((r)) which, 
with (t) in units of mcs, should equal Y 2 +T\n(A 2 ). 
These are for L = 24 for various \H\ values. The 




J/T 



Fig. 1. Mean lifetime, (t), in units of mcs as a function of 
T _1 at \H\=3J/A for the square lattice using the Glauber 
dynamic. The symbols are s=3 MCAMC data. The dashed 
line is the low-temperature prediction with the prefactor 
set to one, while the heavy solid line is with the correct 
prefactor of 3/16 [16]. The vertical dashed line marks T c . 
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Fig. 2. Average lifetime, Tln((r}), in units of J and (t) 
in units of mcs, as a function of T for a L = 24 square 
lattice using the Glauber dynamic. The symbols at finite T 
are s=3 MCAMC data. The values of \H\/J used, reading 
from the top, are 5/6, 0.90, 0.95, 1.00, 1.05, and 1.10. The 
symbols at T = are the exact values of Y2 , and the lines 
use the exact prefactor A2. The prefactor is 3/8 for li = 2 
(J < |H| < 2J) and 3/16 for £ 2 = 3 (2J/3 < \H\ < J). 
The prefactor at \H\ = 1 is not known. 

intercept is the exact low-temperature predic- 
tion, r 2 , and the linear slope is the exact low- 
temperature prefactor, A 2 . These data indicate 
that both the value of T 2 and the prefactor agree 
at low enough temperatures with the simulation 
data. However, exactly how low is low enough to 
agree with the data depends on the value of \H\. 
In particular, for values far from the special value 



4 




|H|/J 

Fig. 3. The average lifetime, (t), in units of Monte Carlo 
steps (mcs) as a function of \H\ at T = 0.2T C for the simple 
cubic lattice using the Glauber dynamic. The symbols arc 
both MCAMC data using up to 3 transient states and 
projective dynamics simulations [19] . The heavy solid line 
is the low-temperature prediction with the prefactor set to 
one, while the dashed line segments are with the correct 
prefactor [16]. 

where £i changes, the MCAMC values agree with 
the low-temperature predictions even for tempera- 
tures as large as T = J/2 (note T c « 2.26J). How- 
ever, closer to the special value where £2 changes 
much lower temperatures are needed before agree- 
ment with the low-temperature predictions are 
seen. Hence the discontinuity in (r) due to the dis- 
continuous prefactor is not seen in the simulation 
data. 



strong fields. No prefactors for these values (since 
\H\ > 2 J) are predicted by [16]. For 4 J < \H\ < 
6 J the prefactor can be evaluated using absorbing 
Markov chains, and is 7/6. This result is consistent 
with the slope for T < 0.2 J. 

The most striking result shown in Fig. 4 is that 
for both \H\ = 3.5J and \H\ = 3.9 J the MCAMC 
results do not agree with the low-temperature pre- 
dictions [16-18] . These are shown by the lowest cor- 
responding symbols at T = 0. An analysis of the 
path taken by the nucleating droplet shows that 
there is a higher saddle point than that predicted 
by the low-temperature results of Eq. 5. The value 
expected at this higher saddle point is shown by 
the higher corresponding symbols at T — 0. For 
3J < \H\ < 4J this saddle has energy 28J - 6\H\, 
corresponding to an L-shapcd droplet with three 
overturned spins. The region of \H\ for which there 
is a higher saddle point shrinks with increasing £ 3 , 
but is always finite near any non-zero value of \H\. 
This higher saddle point removes the discontinu- 
ity in T 3 near the values of \H\ where £3 changes. 
In fact, T 3 is continuous for the corrected saddle, 
just as T 2 is continuous. The discontinuity in (r) 
due to a discontinous prefactor has not yet been 
supported by the MCAMC data for the sc lattice. 
One might expect the same behavior as seen in the 
square lattice, which keeps the discontinuity from 
being seen at any finite temperature. 



5. Cubic-lattice Results 

Figure 3 shows results for the sc lattice. The 
MCAMC and projective dynamics [19] results both 
agree with each other, and away from the values 
of I H I where £3 change they agree with the low- 
temperature prediction with the predicted prefac- 
tor [16]. At these special \H\ values T3 decreases 
by 4 J. This would predict that near these values as 
\H\ decreases the lifetime of the metastable state 
would decrease. This is not physically intuitive, 
and the simulation data do show evidence of the 
predicted discontinuity for T = 0.2T C . 

Figure 4 shows results of MCAMC simulations of 
the sc Ising model at low temperatures for various 



6. Conclusions 

The Monte Carlo with Absorbing Markov 
Chains (MCAMC) algorithm for long-time simu- 
lations of dynamics was described. The descrip- 
tion was for the sc Ising ferromagnet, but can be 
generalized to other discrete systems. The s = 1 
MCAMC algorithm corresponds to the n-fold way 
algorithm [4] in discrete time [5]. In many cases 
the s — 1 MCAMC can require many orders of 
magnitude less computer time than conventional 
dynamic simulations to perform the same calcula- 
tion. The s > 1 MCAMC can require many orders 
of magnitude less computer time than the s = 1 
MCAMC. These exponential decreases in com- 
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Fig. 4. For the simple cubic lattice with Glauber dynamics, 
Tin ((t)) is shown for several values of \H\. The symbols at 
finite T are MCAMC data using up to s = 4. For \H\ = 2.7 J 
below about T = 0.7 J the MCAMC results seem to agree 
with a linear approach to the exact prediction, indicated 
by the same symbol at T = 0. The same is also true for 
the results at \H\ = 5 J where the T = intercept is equal 
to 2. However, for both \H\ = 3.5J and \H\ = 3.9J the 
low-temperature predictions of Eq. 5 are not valid and do 
not fit the MCAMC results. These invalid results are shown 
by the lower corresponding symbols at T = 0. Rather, 
the MCAMC results tend toward the corrected predictions 
given by the higher corresponding T = symbols. 



puter time are accomplished without changing the 
underlying dynamic. 

The MCAMC algorithm has been applied 
to check the low-temperature predictions for 
metastablc decay in both the square and sc Ising 
ferromagnet. In particular, the predicted discon- 
tinuity in the average lifetime at the values of 
| if | where the size of the droplet at the saddle 
point changes discontinuously was investigated. 
It was found that for the square lattice, where 
the discontinuity is only in the prefactor, a finite 
temperature simulation sees no evidence of this 
prefactor discontinuity. Instead, the temperature 
at which the simulations must be performed to 
see the low-temperature predictions decreases as 
the discontinuity is approached. For the sc lattice 
a discontinuity in both the prefactor and in the 
energy of the nucleating droplet was predicted 
[16-18]. No evidence for these discontinuities were 
found in the simulation data. The 'exact formula' 
for the cubic lattice nucleating droplet, Eq. 5, was 
rather found to be incorrect near the discontinu- 
ities. The corrected exact formula has no disconti- 



nuity in the energy of the nucleating droplet. 
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